DESY 10-037, RM3-TH/10-6, ROME1/1468, ROM2F/2010/06 



Non-perturbative renormalization of quark bilinear 
operators with Nj = 2 (tmQCD) Wilson fermions 
and the tree-level improved gauge action 




M. Constantinou^, P. Dimopoulos^), R. Frezzotti( c ' d \ G. Herdoiza^ 
K. Jansen^ V. Lubicz^, H. Panagopoulos^, G.C. Rossi( cd \ 
S. Simula^), F. Stylianou^ , A. Vladikas^ 

W Department of Physics, University of Cyprus, P.O.Box 20537, Nicosia CY-1678, Cyprus 

( b ) Dip. di Fisica, Universita di Roma "La Sapienza" , and INFN, Sezione di Roma, 
P.lc A. Moro 2, 1-00185 Rome, Italy 

( c ) Dip. di Fisica, Universita di Roma "Tor Vergata" , 
Via della Ricerca Scientifica 1, F00133 Rome, Italy 

( d ) INFN, Sez. di Roma "Tor Vergata", c/o Dip. di Fisica, Universita di Roma "Tor Vergata", 
Via della Ricerca Scientifica 1, 1-00133 Rome, Italy 

^ NIC, DESY, Platanenallee 6, D-15738 Zeuthen, Germany 

W Dip. di Fisica, Universita di Roma Tre, Via della Vasca Navale 84, 1-00146 Rome, Italy 

(a) INFN, Sez. di Roma III, Via della Vasca Navale 84, 1-00146 Rome, Italy 

Abstract 

We present results for the renormalization constants of bilinear quark operators 
obtained by using the tree-level Symanzik improved gauge action and the Nf = 2 
twisted mass fermion action at maximal twist, which guarantees automatic 0(a)- 
improvement. Our results are also relevant for the corresponding standard (un- 
twisted) Wilson fermionic action since the two actions only differ, in the massless 
limit, by a chiral rotation of the quark fields. The scale-independent renormalization 
constants Zy, Z& and the ratio ZpjZg have been computed using the RI-MOM 
approach, as well as other alternative methods. For Za and Zp/Zs, the latter are 
based on both standard twisted mass and Osterwalder-Seiler fermions, while for Zy 
a Ward Identity has been used. The quark field renormalization constant Z q and 
the scale dependent renormalization constants Z$, Zp and Zp are determined in the 
RI-MOM scheme. Leading discretization effects of 0(g 2 a 2 ), evaluated in one-loop 
perturbation theory, are explicitly subtracted from the RI-MOM estimates. 



1 Introduction 



Non-perturbative renormalization is an essential ingredient of lattice QCD calculations 
which aim at a percent level accuracy. In this paper, we present calculation methods and 
results for the renormalization constants (RCs) of bilinear quark operators for the lattice 
action used by the ETM Collaboration. The simulations with Nf = 2 dynamical flavours 
employ the tree-level Symanzik improved gauge action and the twisted mass fermionic 
action at maximal twist. In the chiral limit, the latter is related to the standard Wilson 
fermionic action by a chiral transformation, under which quark composite operators behave 
in a definite and simple way. Therefore the RCs we compute here can be also employed to 
renormalize bilinear quark operator matrix elements computed with the same glue action 
but standard (un-twisted) Wilson quarks. 

For the computation of the bilinear operator RCs we have used the RI-MOM method 
PP. For the calculation of the scale independent RCs, Zy, Za and the ratio Zp/Zs, a set of 
alternative methods have also been used. In particular, we have determined Zy by imposing 
the validity of the axial Ward identity whereas for Za and ZpjZ$ we have employed a new 
method which makes a combined use of standard twisted mass (tm) and Osterwalder-Seiler 
(OS) formalism. A comparison of the results obtained with different approaches provides 
an estimate of the size of the systematic errors affecting the calculation of the RCs at fixed 
lattice bare coupling. Needless to say, in renormalized quantities the discretization errors 
coming from both RCs and bare quantities will be removed by continuum extrapolation, 
as illustrated with an example in sect. 3.3. 

Throughout the paper we follow the standard practice of labelling the RCs according 
to the notation adopted in the un-twisted Wilson case. Thus, for instance, Zy and Za 
denote the RCs of the local vector and axial- vector currents of the standard Wilson action 
even though, in the maximal twist case, they renormalize in the physical basis the local 
(charged) axial- vector and polar- vector currents respectively. We refer to Table [1] for the 
explicit presentation of the renormalization pattern. 

We have computed the RCs at three values of the gauge coupling, namely (3 = 3.80, 
3.90 and 4.05, which correspond to inverse lattice spacing a -1 ~ 2.0, 2.3 and 2.9 GeV. 
Preliminary results of this work have been presented in ref. [2]. Further details concerning 
the ETMC simulations can be found in refs. [3]- [6]. 

The paper's contents are the following. In Section [2] we present the novel approach of 
computing the scale independent RCs Za and Z P jZg and the Ward identity determination 
of Zy. In Section [3] we discuss the RI-MOM approach applied to tm quarks and provide the 
details of the numerical analysis. A collection of our final results for the RCs can be found 
in Tables HI [5] and [6j In Tables 0] and [5] the predictions of one loop boosted perturbation 
theory [TJ [8] for the various RCs are also given. Finally in the Appendix we present the 
proof of the automatic (9(a)-improvement of the RI-MOM RCs calculated with maximally 
twisted quarks. 

An RI-MOM calculation of the RCs of bilinear quark operators similar to the one de- 
scribed in the present paper but for the standard Wilson's plaquette gauge action and non- 
perturbatively improved Wilson (clover) fermions has been recently presented in ref. [9]. 
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2 A novel approach to the calculation of the scale 
independent RCs 

In this section we present a calculation of the scale independent RCs, namely Zy, Za and 
ZpjZg. The evaluation of Zy is based on the PCAC Ward identity, which in the quenched 
approximation has led to very precise results (see e.g. refs. [TQl [EE])- On the other hand, 
Za and Zp/Zg are computed from two-point correlators only with a new method, which 
is based on the simultaneous use of two regularizations of the valence quark action. One 
is the standard tm action at maximal twist, while the other is the OS variant [12]. In the 
so called physical basis these actions can be written in the form: 

Svai = a A J2 J2 Q( x ) (7V - «75 r q W cr + ^ q(x) , (1) 

x q=u,d 

with W cr = — I J2/x V* V M + M cr . The Wilson parameters are r u = — = 1 for the standard 
tm case and r u = = 1 for OS quarks. In the sea sector we have two degenerate quarks, 
regularized in the standard tm framework. 

Let us illustrate the general idea of the calculation of the scale independent RCs, Za 
and Zp/Zs- It is based on the observation that the axial transformations of the quark 
fields from the physical basis to the so called twisted basis (primed fields), 

q = exp (i^r q ix/A) q' and q = q' exp (ij 5 r g 7r/4:) , (2) 

transform the tm and OS actions of Eq. (pQ) into an action with the Wilson term in the 
standard form (i.e. having no 75 and r u = — 1) and a mass term which takes the form 
(ifi u u >n f5u' — ifi^d'^d'j in the tm case and {i^i u u'^u' + i^dd'^d'^j in the OS case. This also 
implies, in turn, that the RCs for the operators in the twisted basis, being defined in the 
chiral limit, are the same for the Wilson, tm and OS cases. Consider, now, a non-singlet 
quark bilinear operator, Or = uTd, defined in terms of the fields of the physical basis. 
Under the axial transformations of Eq. (j2j) this operator transforms into an operator in 
the twisted basis, denoted as Op and O r for the tm and OS case respectively. In general, 
the two operators are not of the same form. However their renormalized matrix elements 
between given physical states, say Z Qt (a|O r \(3) tm and Zo {a\Of\j3) os are estimates of the 
same physical matrix element (a|(Or)i?|/3) up to 0(a 2 ) errors: 

(a\(O r ) R \P) = Z 0p (a\Or\(3) tm + 0(a 2 ) = Z Qf («|O r |/3)° 5 + 0(a 2 ) . (3) 

RCs are named, as anticipated, after the twisted basis, in which the Wilson term has its 
standard form. The operator renormalization pattern in the physical and twisted bases is 
given in Table [1] for both OS and tm formulations at maximal twist. The primed operators 
refer to the twisted basis while the unprimed ones to the physical basis. Notice that the 
condition of maximal twist ensures that cut-off effects are of order 0(a 2 ) [T5] . 

The main point of Eq. ([3]) is that the tm and OS determinations of the same continuum 
matrix element are equal up to discretization effects. Our proposal for the determination 
of Zp/Zs and Za, described in the following subsections, is based on this observation. 
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ZrT^v^d 


= iZ T T^ ud 


(Tr) 


[iv,ud ~ 


- Zj-T^y^d ~ 


= ZrP^^d 



Table 1: Renormalization pattern of the bilinear quark operators for the OS and tm case at 
maximal twist. The unprimed operators refer to the physical basis while the primed ones 
to the twisted basis. The symbols P u d and T^^d indicate the operators u^d and ua^^d. 

2.1 Calculation of Zp/Z$ 

The method is based on comparing the amplitude g w = (0\P\n), computed both in tm and 
OS regularizations. We start by considering, in the physical basis, the correlation function 

C PP {t) = - J2(^d(x) d l5 u(0)) , (4) 

X 

which at large time separations behaves like 

C PP (t) ~ ^£[exp(-rM) + exp(-m w (T - t))) (5) 

From Tabled] we see that in the twisted basis this corresponds to C s > S '(t) (for OS regular- 
ization) and C P > P <(t) (for tm regularization) : 

[C PP (t)] cont = Z P [Cp, p/ (t)] tm + 0(a 2 ) = Z 2 s [C S 's'(t)] os + 0(a 2 ) , (6) 

which implies 

[g n r nt = Zp[g n ] tm + 0(a 2 ) = Z s [g n ] os + 0{a 2 ) . (7) 

The ratio Zp/Zg is extracted from the above equation in the chiral limit. As our simulations 
are performed at finite quark masses, an extrapolation of our Zp/Z^-estimators to zero 
quark mass is necessary. 

2.2 Calculation of Za 

In order to compute Za-, we consider the calculation of the charged pseudoscalar meson 
decay constant /„., in both OS and tm regularizations. We start with the tm regularization. 
From the axial Ward identity in the physical basis, we have the well-known result for the 
pseudoscalar meson decay constant [Tl"] : 



[ur nt = [ur + o(a 2 



o 1 t m 

ml 



+ 0(a 2 



Note that in this case no RC is needed [H]. Thus the pion decay constant can be extracted 
from the large time asymptotic behaviour of [Cprpi(t)] tm . 
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In the OS formulation we consider, besides the correlator Cpp of Eq. (j3J), also Cap, 
which is defined by 

C AP (t) = ]>>7o7 5 d(x) d 75 u(0)) . (9) 

x 

Its large time asymptotic behaviour is: 

C AP (t) ~ ^[exp(-m^) - exp(-m^(T - £))] , (10) 

with 

Up = (0|Ao|7r)(7r|P|0) = [f T m n tg n ] os . (11) 

The last equation may be solved for /„., which is thus obtained from the correlation func- 
tions of Eqs. (J5J and (ITUj) . With the aid of Table [TJ this solution is expressed in terms of 
OS-regularized quantities in the twisted basis (note Up = ^U'S') as follows: 

[/.r 4 = Z A \f K )°* + 0(a*) = Z A [g J? { ^ ]OS +0(a>) . (12) 

Combining Eqs. (JH]) and ( |T2l) . we obtain 

[A] con * = iU} tm + 0(a 2 ) = Z A [U] os + 0(a 2 ) , (13) 

from which an estimate of Z A can be extracted. Again, results obtained at finite quark 
masses need to be extrapolated to the chiral limit. In ref. [15] the discretization effects 
affecting the quenched pseudoscalar decay constant in the tm and OS regularizations have 
been studied in detail. The tm and OS results were compatible within one standard 
deviation for three values of the lattice spacing, indicating the smallness of 0(a 2 ) cut-off 
effects. 



2.3 Calculation of Zy 

The determination of the renormalization constant Zy can be done using solely tm quarks. 
It is based on the comparison of the point-like axial current A^^, which renormalizes with 
Zy (see Table [TJ, with the exactly conserved one-point split current [V^]^ s = [^]^ S - 
The four-divergence of the latter is exactly equal to the sum of the valence quark mass 
values, (fii + fi-i), times the corresponding pseudoscalar density. Therefore, we extract Zy 
by solving the equation: 

Zyd [Cy, P ,(t)] tm = (/i! +/i 2 ) [Cp, P ,(t)} tm (14) 

where do is the symmetric lattice time derivative, and extrapolating results to the chiral 
limit. 



2.4 Results 

In Table [2] we give various details of our Nf = 2, partially quenched simulations. The 
smallest sea quark mass corresponds to a pion of about 300 MeV for the two larger /3 val- 
ues, while for the smallest coupling the lightest pion weights approximately 400 MeV. The 
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L 3 x T 






N f 

1 v con j 


A 2 


3.80 


24 3 x 48 


0.0080 


{0.0080, 0.0110, 0.0165, 0.0200} 


400/240 


A 3 






0.0110 


{0.0080, 0.0110, 0.0165, 0.0200} 


400/240 


A, 






0.0165 


{0.0080, 0.0110, 0.0165, 0.0200} 


400/240 


B 1 


3.90 


24 3 x 48 


0.0040 


{0.0040, 0.0064, 0.0085, 0.0100, 0.0150} 


240/240 


B 2 






0.0064 


{0.0040, 0.0064, 0.0085, 0.0100, 0.0150} 


240/240 


B 3 






0.0085 


{0.0040, 0.0064, 0.0085, 0.0100, 0.0150} 


240/240 


B 4 






0.0100 


{0.0040, 0.0064, 0.0085, 0.0100, 0.0150} 


240/240 


Ci 


4.05 


32 3 x 64 


0.0030 


{0.0030, 0.0060, 0.0080, 0.0120} 


130/240 


c 2 






0.0060 


{0.0030, 0.0060, 0.0080, 0.0120} 


130/160 


c 3 






0.0080 


{0.0030, 0.0060, 0.0080, 0.0120} 


130/160 



Table 2: Details of the simulations performed for computing RCs. The number N con f of 
gauge configurations we analysed is given for both the case of the scale independent RCs 
discussed in this section (first figure) and the RI-MOM analysis of Sect. 3 (second figure). 

highest sea quark mass is around half the strange quark mass. For the inversions in the 
valence sector we have made use of the stochastic method (one-end trick of ref. [16]) in or- 
der to increase the signal to noise ratio. Propagator sources have been placed at randomly 
located timeslices. This turned out to be an optimal way to further reduce the autocorre- 
lation time. Our correlation functions are computed for all combinations of valence quark 
masses appearing in Table |2J In order to deal with effectively independent measurements 
of RC-estimators, we have selected gauge field configurations separated by 20 HMC (length 
1/2) trajectories. Within each ensemble of such gauge configurations statistical errors are 
evaluated using a jackknife procedure. For results involving an extrapolation in the sea 
quark mass and/or combined fits at several lattice couplings, statistical errors have been 
estimated with a bootstrap procedure for a 1000 bootstrap events. Typical plots on the 
data quality of Zy, Za and Zp/Zs are shown in Fig. [TJ in this example we show the case 
a f i Tea = = a / i 2 for the three values of the gauge coupling. 

Various methods have been implemented in taking the chiral limit (in the sea and 
valence quark sector). A first method consists in calculating the RCs at fixed sea quark 
mass for three choices of valence quark pairs, namely with a/ii = afi 2 , or for all pairs of 
(a/xi,a/i2) or for (au.i,afi 2 ) > afi sea and taking the "valence chiral limit" using a linear fit 
in terms of the sum of the valence quark masses. Subsequently, the RCs were quadratically 
extrapolated to the sea quark chiral limit. Note that a quadratic dependence on a\i sea is 
expected from the form of the sea quark determinant, assuming that lattice artefacts on the 
RCs are not sensitive to spontaneous chiral symmetry breaking. However we have verified 
that a linear fit in fi sea leads to compatible results, albeit with larger final errors. A second 
method consists in inverting the order of the two chiral limits. A third method is simply 
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(c) 

Figure 1: The plateaux quality for Z v (a), Za (b) and ZpjZg (c) for three gauge couplings 
and masses a/i™ n = a//i = a/i 2 . From Table [2| we have au,™™{fi = 3.80) = 0.0080, 
= 3.90) = 0.0040 and a//™ n (/3 = 4.05) = 0.0030. iVote tiat in Fig.(b) the data 
have been shifted for clarity. 



the extraction of the RCs from the subset of data satisfying /J, va i en ce = Usea- This allows 
reaching the chiral limit with a single, linear extrapolation in the quark mass. In most 
cases the quality of the fits is very good and the final results, obtained from these different 
extrapolations, are compatible within one standard deviation. We have opted to quote 
as final results those produced by the first method, for all pairs of (afii, 0,0-2), followed by 
a linear fit in terms of a 2 o, 2 sea . In Fig. [2] (left panel) we show the "valence chiral limit" 
extrapolation of the three scale independent RCs, computed at the lightest sea quark mass 
for each gauge coupling. In Fig. [5] (right panel) we present the chiral limit in the sea sector 
for all the RCs for each gauge coupling. 

Our final results are collected in Table H] (see rows labelled as "Alt. methods"). For 
comparison with the RI-MOM results and conclusions on the precision of the methods 
used, see below. Here we point out that the Zy value at (3 = 3.90 compares nicely with 
the result of Ref. [17], produced from 3- and 2-point correlation functions, in the context 
of a calculation of the pion form factor. 
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Figure 2: Left panel: valence chiral limit for the three scale independent RCs at (3 = 3.80, 
a/x sea = 0.0080 (a), (5 = 3.90, a/i sea = 0.0040 (c) and (3 = 4.05, a/z sea = 0.0030 (e). Right 
panel: sea chiral limit for the same constants at (3 = 3.80 (b), f3 = 3.90 (d) and (3 = 4.05 
(f). Note that the data for Zp/Zs have been shifted for clarity. 
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3 Renormalization constants in the RI-MOM scheme 



The non-perturbative determination of the RCs in the RI-MOM scheme pQ is based on 
the numerical evaluation, in momentum space, of correlation functions of the relevant 
operators between external quark and/or gluon states. In this work we are interested in 
the case of the bilinear quark operators Or = uTd, where T = S,P, V, A, T stands for 
1, 75; T/ttj TmT5> respectively. The fields u and d of the operator Or belong to a twisted 
doublet of quarks, regularized by the tm action of Eq. (JTJ, with r u = —r^ = 1. 
The relevant Green functions are the quark propagator 

S q (p) = a 4 £ e"** (q(x)q(O)) (? = «, d) , (15) 

X 

the forward vertex function 

Grip) = a 8 £ e-*-^ («(x)O r (0) %)) , (16) 



and the amputated Green function^ 

Ar(p) = SM^GrWSiip)- 1 . (17) 

The RI-MOM renormalization scheme consists in imposing that the forward amputated 
Green function, computed in the chiral limit in the Landau gauge and at a given (large 
Euclidean) scale p 2 = /i 2 , is equal to its tree- level value. In practice, the renormalization 
condition is implemented by requiring in the chiral limit that 

Z- x Z T V f {p)\^ = Z-%Tr[A f (p)P f ]| p2=Al2 = 1 , f = e -^/4 re *W4 > (is) 

where Pf is a Dirac projector satisfying Tr [f Pf] = lj^l The quark field RC Z q , which 
enters Eq. (fT8l) . is obtained by imposing, in the chiral limit, the conditional 



Z' 1 S 1 (p)| p2=M 2 = Z' 1 — Tr 



p 2 



1. (19) 



In order to minimize discretization effects, we select momenta with components p u = 
[2n/ L u ) n u in the intervals 

•H KmImIm! ' fMi = 24 (/! = 3 ' 8a,ld3J) (20) 



^^The calculation of Gr(p) in Eq. (|16[) involves the propagator 5^(0, y) which, in the tm approach, equals 
7 5 S' u (y,0) t 7 5 . Similarly, in Eq. ([17]) the propagator S d (p) stands for J2 y e w ' v S d (0,y) = 7 5 S'„(p) t 7 5 . 

2 The appearance of T in Eq. (|18l) is a trivial consequence of the fact that RCs, as we said, are named 



after the (here unphysical) twisted quark basis while the operator Or = uTd is expressed in the physical 
quark basis. 

3 Strictly speaking, the renormalization condition of Eq. ([T^| defines the so called RI' scheme. In the 
original RI-MOM scheme the quark field renormalization condition reads 



V^Tr 



7m" 



9pn 



1 . 



The two schemes differ in the Landau gauge at the N 2 LO. In this work, when perturbative predictions are 
used, the difference between the two schemes has been properly taken into account. 
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and 



• ( ^« 5) (21) 



(L y is the lattice size in the direction v). The time component of the four- momentum 
is shifted by the constant Ap^ = ir/L^, in order to account for the use of anti-periodic 
boundary conditions on the quark fields in the time direction. Furthermore, we have only 
considered for the final RI-MOM analysis the momenta satisfying 



A 4 (p) = \ 2 < 0.28, (22) 



where 



p u = - sin(a p u ) , (23) 
a 

which helps in minimizing the contribution of Lorentz non-invariant discretization effects 
(cfr Eq. d3U). 

In order to improve the statistical accuracy of the RCs, we have averaged the results 
obtained from the correlation functions of the operators Or = uTd and 0' T = dTu. Sim- 
ilarly, in the case of the quark field RC, we have computed Z q by averaging the results 
obtained for the up and down quark propagators. 

The RCs, calculated in the chiral limit in the way described above, are automatically 
improved at 0(a). Actually in the maximal twist situation 0(a) improvement holds for 
all external momenta p and non-zero (valence and sea) quark mass at the level of the basic 
quantities Tr[A r (p) P r ] and Tr [(J> S g (p)~ 1 ) /p 2 ] entering in Eqs. (!TB ]) - (fl9"]) . A proof of this 
statement, which follows from an analysis based on the symmetries of the tm action at 
maximal twist and the 0(4) symmetry of the underling continuum theory, is given in the 
appendix^ 

Details of the RI-MOM simulation are collected in Table |2j The lattice parameters are 
those used also in the determination of the scale independent constants with the methods 
discussed in the previous section. However, a different ensemble of independent gauge 
configurations has been analyzed in each case. Moreover, the inversions in the valence 
sector for the RI-MOM study have been performed without using the stochastic approach, 
but with point-like sources randomly located on the lattice for each gauge configuration. 



3.1 Analysis of the twisted mass quark propagator 

The necessary Green functions for the RI-MOM determination of RCs of bilinear quark 
operators are the quark propagator S q {p) and the amputated vertex Ar(p), evaluated in 
momentum space in the Landau gauge. In this section, we first illustrate the results on 
the tm quark propagator. 

At 0(a), the spin-flavour structure of the lattice artefacts of the quark propagator 
differs from that of the standard Wilson case, due to the twisted Wilson term in the 

4 With the standard Wilson or Clover actions, the RCs determined with the RI-MOM method are 
also 0(a)-improved. The RC-estimators are however improved only at large momenta, where spontaneous 
chiral symmetry breaking effects can be neglected [15] . This is sufficient, as the RI-MOM scheme is defined 
precisely in the high-momentum region. 



9 



action. With tm fermions, the explicit breaking of parity at finite lattice spacing allows 
for the presence of a parity violating contribution in the quark propagator. By neglecting 
0(4) symmetry violating effects, which only appear at C(a 2 ), the inverse quark propagator 
can be expressed in terms of three scalar form factors, as follows: 

S.ip)- 1 = i^p 2 ) + S 2 (p 2 ) -i 15 Mp 2 ) ■ (24) 

At large p 2 , Si and E 2 are related to the quark field RC Z q (see Eq. ( 1T9|) ). and to the 
renormalized quark mass fi q respectively [19]. The parity violating term proportional to 
E 3 represents an 0(a) discretization effect, induced by the twisted Wilson term in the 
actionjf] At maximal twist, the form factors S 1>23 are given at tree-level by 

S 1 (p 2 ) = l , MP 2 )=N , £ 3 (p 2 ) = ^V. (25) 

In Fig. [21 the non-perturbatively determined form factors S 1)2j3 are plotted against a 2 p 2 , 
for (3 = 3.90 and a[i sea = 0.0040. The following observations are in order: 

• The form factor Si exhibits a nice plateau at large p 2 , with a tiny slope which is 
partly due to a NLO renormalization scale dependence (the quark field anomalous 
dimension vanishes at LO in the Landau gauge) and partly to 0(a 2 ) discretization 
effects. The quark field RC Z q is obtained from the Si results at large p 2 , extrapolated 
to the chiral limit (see Eq. ffT9"j) ). 

• The form factor S 2 also exhibits a plateau at large p 2 , which is expected to be 
proportional, up to the quark field RC, to the renormalized valence quark mass fi q 
in the RI-MOM scheme. Indeed, at variance with the cases of S x and S 3 , a clear 
separation among the results obtained for S 2 at different bare valence quark masses 
is visible in the plot. 

• As expected, the form factor S3, which represents a pure 0(a) discretization effect, 
has opposite sign for the up and down quark propagators {r u = —r^ = 1). Its 
behaviour is quite close to the tree level estimate of Eq. (|25|) . We find S 3 ~ ± c q ap 2 /2, 
with c q ~ 0.8 ^ 0.9. 

3.2 Renormalization constants 

The RCs in the RI-MOM scheme are determined by closely following the procedure sum- 
marised in section [3] above and described in detail in ref . [18] . This procedure involves two 
main steps: the chiral limit extrapolation of the RCs, at fixed coupling and renormalization 
scale, and the study of the renormalization scale dependence. 

3.2.1 Chiral extrapolations 

RI-MOM is a mass independent renormalization scheme. Since in practice the RCs are 
obtained at non vanishing values of the quark masses, an extrapolation of the results to 
the chiral limit, both in the valence and the sea quark masses, must be performed. 

5 In the standard Wilson case, the analogous 0(a) artifact contributes to the form factor £2- 
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Figure 3: The form factors Si (top), S 2 (center) and S 3 (bottom), computed at (3 = 3.90 
with a/z sea = 0.0040, as a function of a 2 p 2 . 
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The validity of the RI-MOM approach relies on the fact that, at large p 2 , Green func- 
tions are expected to depend smoothly on the quark masses, since their non-perturbative 
contributions vanish asymptotically in that limit. However, specific care must be taken 
in the study of the pseudoscalar Green function Vp, since the leading l/p 2 -suppressed 
contribution is divergent in the chiral limit, due to the coupling with the Goldstone bo- 
son [TJ [201 EI] • Moreover, for tm fermions, the explicit breaking of parity at finite lattice 
spacing induces a coupling of the Goldstone boson also with the scalar operator. Though 
the latter is suppressed at 0(a 2 ), its effect may be not negligible at the couplings considered 
in the present simulations. 

In order to subtract the pseudoscalar mass dependence of the amputated vertex func- 
tions Vp and Vs, at each given p 2 we fit their values at different quark masses to the 
Ansatz 



V P(5 )(p 2 ,/ii,/i 2 ) = A P(s) (p 2 ) + B P{s) (p 2 ) M* g fa,to) + A Z P{ p iPj , , (26) 



where Mps(/zi,/z 2 ) is the mass of the pseudoscalar meson composed by valence quarks of 
mass iii and /i 2 (with r x = — r 2 ) %. The first term in the fit, i.e. the function A P ( S )(p 2 ), 
represents the Green function Vprs) i n the (valence) chiral limit, from which the RC Zp(g) is 
eventually extracted. The last term in Eq. ( 126]) accounts for the presence of the Goldstone 
pole. Using the results of the fit, we can define the subtracted correlators 

W, ,*,,*> = V P{S ^,^) - sggg) • (27) 

The effect of the subtractions is illustrated in Fig. |U We find that, while in the case of 
Vp the contribution of the Goldstone pole is clearly visible, the subtraction is practically 
irrelevant in the case of Vs, indicating a strong 0(a 2 ) suppression of the parity violating 
coupling. We have also verified that, the alternative procedure [22] for the Goldstone pole 
subtraction, based on the definition 

,,sub / 2 \ M PS (m,II 2 ) Vp (5 )(p 2 ,/ii,/i 2 ) - M PS {fi 3 ,fI 4 ) Vp( 5 )(p 2 ,/i3,^4) 

Vp(S)\P > ^1^2, ^3,^4 J — 772—7 7 772—7 7 • 

1 ; Mp s {ni, fi 2 ) - M PS (fi 3 , /i 4 ) 

(28) 

leads to consistent results in the chiral limit, within our statistical errors. 

The dependence of Z v , Za and Z T , which are not affected by the Goldstone pole 
contamination, on the valence quark masses, is shown in Fig. [51 for all three couplings. 
For illustration purposes, a specific value of the sea quark mass has been chosen in each 
case. We see that the valence quark mass dependence of the RCs is rather weak and well 
consistent with a linear behaviour. 

In the tm formulation of lattice QCD with Nf = 2 flavours of degenerate sea quarks, 
the fermionic determinant is a function of the sea quark mass squared. Therefore, we 
expect that the RCs, which by definition are insensitive to the effects of spontaneous chiral 
symmetry breaking, will also exhibit the same dependence. This is shown in Fig. [61 We 



6 In our preliminary analysis presented in [2], the fitting function used to subtract the Goldstone pole 
contribution was expressed in terms of the sum of quark masses [i\ + /12 rather than the pseudoscalar 
meson mass squared Mps(ni, /^2) 2 j as in Eq. (|26p. Since the latter is simply proportional to fa + /i 2 at 
maximal twist, the two fits lead to completely consistent results. 
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find that all RCs depend very weakly on the sea quark mass squared. We thus obtain our 
mass independent RCs by performing the sea quark mass extrapolation to the chiral limit 
linearly in a 2 /i 2 ea . We also checked that chiral extrapolations based on a first or second 
order polynomial fit in a\i sea lead within errors to practically equivalent results. 

3.2.2 Renormalization scale dependence and subtraction of the 0(g 2 a 2 ) effects 

Once the RCs have been extrapolated to the chiral limit, we investigate their dependence 
on the renormalization scale by evolving, at fixed coupling, the RCs to a reference scale 
Ho = 1/a. This is done using 

Z T (a, Ho) = C t (ho, h) z r(a, h) , (29) 

where the evolution function C-p is expressed in terms of the beta function (3(a) and of the 
anomalous dimension of the relevant operator 7r(ct) by 



C r (Ho,H) = exp 



a M 7r(a) , 
— — da 



(30) 



This function is known, in the RI-MOM scheme, at the N 2 LO for Zt [23] and at the N 3 LO 
for Z s and Z P [21]. Since Z v , Za and the ratio Z P /Z S are scale independent, they have 
vanishing anomalous dimensions; thus for these quantities we have Cy(hqi n) = 1- 

In the non-perturbative calculation, the RCs evolved to the reference scale Ho = V° 
maintain a dependence on the renormalization scale p 2 at which they have been initially 
computed. We will keep track of this dependence, which (at large enough p 2 ) is mostly due 
to discretization effects, by denoting these RCs as Zr(l/a; a 2 p 2 ), where the first variable 
indicates the renormalization scale, Ho — l/ a 5 an d the second the dependence on the initial 
momentum. 

In Fig. we show the results for all Zr(l/a; a 2 p 2 ) at /3 — 3.9 as a function of a 2 p 2 
(empty symbols). The residual a 2 p 2 -dependence which is observed in these results in the 
large momentum region is practically linear, suggesting that leading discretization effects 
are 0(a 2 p 2 ). As illustrated in the plots, this dependence is particularly pronounced in the 
case of the pseudoscalar RC Z P . 

In order to reduce the size of discretization errors, we analytically subtract from the 
quark propagator and the amputated vertex functions the 0(g 2 a 2 ) contributions, recently 
computed in lattice perturbation theory [SJ. Up to 1-loop and up to 0(a 2 ), the amputated 
projected Green functions Vr(p), defined in Eq. (j!8p . have the simple and general expression 



Vr(p) pert - = 1 + ^-{b? + bP ln(aV) 



+ a 2 Uc«+cf Ma 2 p 2 )) + cP^}} + 0(a*g 2 , g 4 ) . (31) 

The values of the coefficients and Cp are collected in Table [3] for the lattice action used 
in the present study, namely the tree-level Symanzik improved gluon action in the Landau 
gauge and the tm fermionic action at maximal twist 0. 



7 The same results are also valid for the standard Wilson fermionic action, except for the exchange of 
the values of the vector (V) and axial (A) coefficients. 



16 



r 


6 r i; 




r y ' 
C r 




C r 


V 


-0.48369852(8) 





1.5240798(1) 


-1/3 


-125/288 


A 


3.57961385(3) 





0.6999177(1) 


-1/3 


-125/288 


S 


0.58345905(5) 


-3 


2.3547298(2) 


-1/4 


1/2 


P 


8.7100837(1) 


-3 


0.70640549(6) 


-1/4 


1/2 


T 


0.51501972(6) 


1 


0.9724758(2) 


-13/36 


-161/216 



Table 3: Values of the coefficients and Cp entering the 1-loop expression / f3l]) of the 
amputated projected Green functions Vrip)- Results are presented for the case of the tree- 
level Symanzik improved gluon action in the Landau gauge and the tm fermionic action at 
maximal twist. 



A result similar to Eq. (!3T|) also holds for the form factor Si(p) of the inverse quark 
propagator, from which the quark field RC Z q is evaluated: 



Ei(p) 



pert. 



+ of 



9_ 

12tt 

/(c^ + cf ln(a 2 p 2 )) + cf)^ 



} + 0(a 4 g 2 ,g 4 ) . (32) 



Our definition of this form factor on the lattice, which is equivalent to Eq. ( TT9l up to terms 
ofC(a 2 ), is 



Ei(p) = 



12N(p) 



Tr 



E'^Cp)- 1 )/^ 



(33) 



where the sum J2 P only runs over the Lorentz indices for which p p is different from zero 
and N(p) = J2'„ 1- With this definition, the coefficients of Eq. fl32|) take the values 



&W = -13.02327272(7) ; bf = ; < 

m . , 2.07733285(2) 
c« = 1.14716212(5) + ^-U 



(3) = _!_ . 
q 240 ' 

(2) 73 157/180 



360 



iV(p) 



(34) 



Using Eqs. (l3"Tj) and (|32|) . we then define the corrected amputated Green functions and 
the corrected form factor Si as 



Vr(p) c 



MP) - 



Si(pr r - = S 1 (p)--|^a 2 



,(3) Ep P p 



p 2 



} 



p 2 (c('> + cf ln(aV)) + <f } • ( 35 ) 



which are free of 0(g 2 a 2 ) effects. Note that the 0(a 2 ) terms depend not only on the 
magnitude, J2 P P% but also on the direction of the momentum, p p , as manifested by the 
presence of J2 P P P - As a consequence, different renormalization prescriptions, involving 
different directions of the renormalization scale p p , treat lattice artifacts differently. In 
the numerical evaluation of the perturbative correction in Eq. (135 p . we replaced g 2 with a 
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Figure 7: The RCs Z r (l/a; a 2 p 2 ) at (3 = 3.9, evaluated at the reference scale fi = 1/a, plot- 
ted against the initial renormalization scale a 2 p 2 . Filled squares (empty circles) are results 
obtained with (without) the subtraction of the 0(g 2 a 2 ) discretization effects computed in 
perturbation theory, see Eq. ( 135]) . The solid lines are linear fits to the data. 
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simple boosted coupling [25], denned as g 2 = g 2 /(P), where the average plaquette (P) is 
computed non-perturbatively. 

The effect of the subtraction ( l35l) is also illustrated in Fig. [7j which shows that the 
a 2 p 2 -dependence of the RCs is significantly reduced by the perturbative correction. By 
fitting the RCs as 

Z r (l/a; a 2 p 2 ) = Z r (l/a) + A r a 2 p 2 (36) 

in the large momentum region a 2 p 2 ^ 1 (with Ar just constrained to depend smoothly on g 2 
- see Eq. fl38|) ). we find that the slopes Ar are reduced, after the perturbative subtraction, 
at the level of 10 -2 or smaller for Zy, Za, Zt and Z q . For Z$ we find that the slope, 
which is also about 10 -2 before implementing the correction, increases slightly after the 
subtraction. For Zp, on the other hand, the slope is rather large before the subtraction; 
i.e. Ap ~ 5 • 10~ 2 . As shown in Fig. [71 the effect of the subtraction is beneficial also in this 
case, but inadequate for correcting the bulk of the observed a 2 p 2 -dependence. This is not 
completely unexpected: when discretization effects are large, the subtraction of only the 
leading 0(g 2 a 2 ) terms may be not sufficient to reduce them to a negligible level. Similar 
results, with approximately equal values of the slopes Ar, are obtained at all three values 
of the lattice spacing. 

A useful way to investigate the size of lattice artifacts, which are left after the pertur- 
bative subtraction, consists in comparing the renormalization scale dependence of the RCs 
as observed at different lattice spacings. Specifically, at each value of the lattice spacing 
we fix two common values of the renormalization scale, fi and //, and compute the step 
scaling functions 

S r pi, n;a) = — -. 37 

Z r (a,n') 

If not for discretization effects, the step scaling functions should be independent of the 
cut-off and equal to the evolution functions Cr(/i, fi') of Eq. fl29|) . 

In fig. [8] we show the results obtained for the step scaling functions of the quark bilinear 
operators and of the quark field RCs, as a function of (a/ro) 2 , for two common pairs renor- 
malization scales, (pi 2 ,// 2 ) = (8.5,6.0) GeV 2 and (/i 2 ,// 2 ) = (8.5,7.25) GeV 2 , which both 
lie in the interval of momenta accessible at all values of the lattice spacing. We see from 
these plots that the dependence of the step scaling functions on the lattice cut-off is tiny 
in most of the cases, being clearly visible only in the case of £p. Moreover, for the scale 
independent RCs Zy and Za, a linear extrapolation of the results to the continuum limit 
leads to values which are perfectly consistent with unity, within the statistical errors. For 
the other RCs, the continuum limit of the step scaling functions leads to non-perturbative 
results which are in good agreement with the perturbative determinations of the evolution 
functions Cr(/x, //), computed at the N 2 LO for Z T and at the N 3 LO for Z s , Z P and Z q . 
These results provide evidence that higher order perturbative contributions to the anoma- 
lous dimensions of the various operator are not relevant for describing the renormalization 
scale dependence of the RCs in the region of momenta explored in the present calculation. 

In order to account for the residual discretization effects in the calculation of the RCs, 
we follow two different approaches: 

- Extrapolation method (Ml): after subtraction of the 0(g 2 a 2 ) terms according to 
Eq. (135|) . we extrapolate the RCs linearly to a 2 p 2 — > 0. Specifically, we fit -Zr(l/a; a 2 p 2 ) 
with Eq. (136 p in the large momentum region, 1.0 ^ a 2 p 2 ^ 2.2. The slopes Ar exhibit 
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Figure 8: Values of the step scaling functions Sr(//, //; a) evaluated for (/i 2 ,// 2 ) = 
(8.5,6.0) GeV 2 (red circles) and (,u 2 ,/i /2 ) = (8.5,7.25) GeV 2 (blue squares) at the three 
values of the lattice coupling. The solid lines and the empty markers show the results of a 
linear extrapolation to the continuum limit. For comparison, the corresponding perturba- 
tive estimates of the evolution functions Cr(/i, //) are also shown in the plots (black circles 
and squares). 



only a very mild dependence on the coupling. We parameterize this dependence by 
performing a simultaneous extrapolation at the 3 values of the coupling and writing 
the slopes as 

ArG? 2 ) = A r ( 5o 2 ) + A^o 2 ) • (g 2 - gl) , (38) 

where g is the coupling corresponding to the intermediate value (3 = 3.9. The 
intercept of the extrapolation as determined from the fit, and indicated with Zr(l/a) 
in Eq. f )36|) . provides our final estimate of the RC with the extrapolation method. 
The fit is illustrated, for all RCs, in fig. EB 

- p 2 -window method (M2): in this approach we do not attempt any additional sub- 
traction of discretization effects, besides the perturbative one. The final estimates of 
RCs are obtained by taking directly the results of Z r (l/a;a 2 p 2 ) at a fixed value of 
p 2 (in physical units). In practice, this is done by fitting the RCs to a constant in 
the small momentum interval p 2 G [8.0,9.5] GeV 2 . The idea behind this approach 
is that, once RCs are combined with bare quantities, so as to construct the physical 
observables of interest, the residual 0(a 2 ) effects, which are present in both RCs and 
bare quantities, will be extrapolated away in the continuum limit. 

The two methods are compared below, in a specific example. This consists in the scaling of 
the pseudoscalar meson mass, composed by two degenerate quarks of fixed (renormalized) 
mass. As expected, the two methods give different results at fixed lattice cut-off, but the 
difference disappears in the continuum limit. 
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Figure 9: The RCs Z r (l/a; a 2 p 2 ), evolved to the reference scale fi = 1/a, as a function of 
the initial renormalization scale a 2 p 2 . The solid lines are linear fits to the data. 
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3.3 Analysis of systematic errors 



Before presenting our final results for the RCs, we list and discuss in some detail the main 
systematic uncertainties. 

- Finite size effects: the RCs, determined in the RI-MOM scheme at large momenta, 
are short distance quantities and, as such, should not be affected by significant finite 
size effects. In order to verify this expectation, an additional RI-MOM calculation 
at /3 = 3.90 has been performed, on a 32 3 x 64 lattice. A comparison of these results 
with those obtained on the smaller 24 3 x 48 lattice is illustrated in Fig. [10J We find 
that, for the momentum range in which results on both volumes are available, the 
differences are smaller than statistical errors. 

- Subtraction of the Goldstone pole and chiral extrapolations: the chiral extrapolation 
of the RCs is rather delicate for Zp (and to a lesser extent for Z$), due to the presence 
of the Goldstone pole which has to be subtracted. The uncertainty introduced by 
the subtraction of the Goldstone pole contribution to Z$ and Zp can be estimated 
by comparing the results obtained following two different subtraction approaches: 
the fit of Eq. ( 1261) or, alternatively, the procedure based on Eq. ( |28l) . We find that 
differences are always negligible for Zg, whereas for Zp they are at the level of our 
statistical errors. 

As illustrated in Figs. [5] and [61 the chiral extrapolations in the valence quark mass 
for the other RCs and in the sea quark mass for all RCs are quite smooth. They are 
also well consistent with the expected leading linear and purely quadratic dependence 
on the valence and sea quark masses respectively. Therefore, we estimate that the 
uncertainties associated with these chiral extrapolations are negligible with respect 
to the statistical errors. 

- Discretization effects: the analysis of these effects has been presented in section [3.2.21 
The leading 0(g 2 a 2 ) discretization errors have been accounted for, using the analyt- 
ical results of ref. [8]. With the method denoted as Ml in section 13.2. 2\ after the 
perturbative subtraction, residual 0(a 2 p 2 ) effects are extrapolated away, with the 
Ansatz (|36|) . These contributions are not subtracted, instead, in the approach de- 
noted as M2. Discretization effects affecting both RCs and the bare matrix elements 
will be corrected for, by extrapolating the physical (renormalized) quantity to the 
continuum limit. 

As an example of the results obtained by implementing the two approaches, we show 
in Fig. [LTJ the continuum extrapolation of the pseudoscalar mass squared, computed 
at a fixed value of the renormalized quark mass /u,p. The renormalization of the latter 
has been achieved using the estimates Ml and M2 for Zp given in Tables 141 and lljpl. 
Although the two Zp determinations differ by visible discretization effects at finite 
lattice spacing, particularly at the two coarsest lattice, the continuum limit results of 
the pseudoscalar mass squared are consistent within the two determinations. This is 
evidence that the discretization errors as well as other possible systematic effects in 
the analysis of the RI-MOM RCs are well under control, at least within our current 
statistical errors (say at the level of one standard deviation and possibly even less). 

8 Note that in tm lattice QCD we have = Zp 1 . 
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Figure 11: Scaling of the pseudoscalar meson mass squared, computed at a fixed value of 
the quark mass. The Ml and M2 determinations for Z^ = Zp 1 lead to compatible results 
in the continuum limit. 

3.4 Final results and comparison with perturbation theory 

The final results for the bilinear quark operators and the quark field RCs, obtained with the 
RI-MOM method, are collected in Tables 1 and E (labelled as RI-MOM (Ml) and RI-MOM 
(M2)). These results are compared with those obtained in Sec. [2]for the scale independent 
Zy, Za and Zp/Z$ (Alt. methods). In addition, we give in Tables H] and [5] the predictions 
of one loop boosted perturbation theory (BPT) [3, [8], obtained with two definitions of the 
boosted coupling: the first is g 2 = g 2 / (P), also used in the previous section, and the second 
is based on the one-loop tadpole-improved expansion of log(P) [25J. 

From Tables H] and O we see that the largest deviations between the central values of the 
RI-MOM (Ml) determinations and the alternative determinations of Section [2] are at the 
level of 3-4% or smaller for Za and Zy. These differences are not always accounted for by 
the errors quoted for each result. When statistically significant, these differences provide 
an estimate of the systematic uncertainties affecting the calculations. As discussed in the 
previous section, we expect these uncertainties to be dominated by 0(a 2 ) discretization 
effects. Nevertheless, we do not include these differences in the final errors. The reason is 
that such estimates of discretization effects, though meaningful at finite lattice spacing, are 
of no practical significance when the continuum limit of a physical quantity is eventually 
computed. 

The comparison of the Zp/Z$ results, obtained with the RI-MOM and the alternative 
method, deserves a further comment. The difference between the two results is about 15% 
for the coarsest lattice and 4% for the finest one. The rather big difference noticed in the 
former case may be attributed not only to discretization effects, but also to an imprecise 
tuning of the twist angle to maximal twist at this coupling. This is expected to influence 
substantially the OS pseudoscalar density. 

As shown in Tables H] and the non-perturbative results also lie in the range of one- 
loop BPT estimates, with the exception of Z$ and, more notably, of Zp. We emphasise 
that the accuracy of Zp is crucial in the determination of quark masses, since the bare 
quark mass, computed with maximally twisted fermions, renormalizes with of Zp ; see for 
example ref. [26] . 
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& 


Method 


Zy 


Z A 


Zp/Zs 




RI-MOM (Ml) 


0.604(07) 


0.746(11) 


0.580(17) 


3.80 


RI-MOM (M2) 


0.623(05) 


0.727(07) 


0.653(07) 




Alt. methods 


0.5816(02) 


0.747(22) 


0.498(22) 




BPT 


0.61-0.70 


0.71-0.77 


0.72-0.80 




RI-MOM (Ml) 


0.624(04) 


0.746(06) 


0.626(13) 


3.90 


RI-MOM (M2) 


0.634(03) 


0.730(03) 


0.669(08) 




Alt. methods 


0.6103(03) 


0.743(18) 


0.579(19) 




BPT 


0.63-0.71 


0.72-0.78 


0.74-0.82 




RI-MOM (Ml) 


0.659(04) 


0.772(06) 


0.686(12) 


4.05 


RI-MOM (M2) 


0.662(03) 


0.758(04) 


0.710(07) 




Alt. methods 


0.6451(03) 


0.746(18) 


0.661(21) 




BPT 


0.65-0.73 


0.74-0.80 


0.76-0.83 



Table 4: Values of Z V ,Z A and Z P /Z S , obtained with the RI-MOM methods Ml and M2 
(see text for details), as well as with the methods described in Sec. (labelled as " Alt. 
method"). The predictions of one loop boosted perturbation theory (BPT) are also shown 
for comparison. 






Method 


Zf(l/a) 




Zf(l/a) 


Zf{l/a) 


3.80 


RI-MOM (Ml) 
RI-MOM (M2) 
BPT 


0.603(15) 
0.668(06) 
0.68-0.75 


0.338(10) 
0.437(04) 
0.49-0.60 


0.743(09) 
0.719(06) 
0.68-0.75 


0.755(07) 
0.755(05) 
0.69-0.76 


3.90 


RI-MOM (Ml) 
RI-MOM (M2) 
BPT 


0.615(09) 
0.669(05) 
0.70-0.76 


0.377(06) 
0.447(05) 
0.52-0.62 


0.743(05) 
0.724(03) 
0.70-0.76 


0.758(04) 
0.755(03) 
0.71-0.77 


4.05 


RI-MOM (Ml) 
RI-MOM (M2) 
BPT 


0.645(06) 
0.678(04) 
0.72-0.78 


0.440(06) 
0.480(04) 
0.55-0.65 


0.763(06) 
0.752(04) 
0.72-0.78 


0.782(05) 
0.778(03) 
0.73-0.79 



Table 5: Values of Z$, Zp and Zj-, obtained with the RI-MOM methods Ml and M2 (see 
text for details). The predictions of one loop boosted perturbation theory (BPT) are also 
shown for comparison. 
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Z v 


Z A 


Z S AS (2 GeV) 


Zf E {2 GeV) 


Z™ s (2 GeV) 


Zf s {2 GeV) 


3.80 


0.604(07) 


0.746(11) 


0.734(18) 


0.411(12) 


0.733(09) 


0.745(07) 


3.90 


0.624(04) 


0.746(06) 


0.713(10) 


0.437(07) 


0.743(05) 


0.751(04) 


4.05 


0.659(03) 


0.772(06) 


0.699(06) 


0.477(06) 


0.777(06) 


0.780(05) 



Table 6: Final results for the RCs of bilinear quark operators and the quark Geld, obtained 
with the RI-MOM method. The scheme dependent RCs are given in the MS scheme at 
scale n=2 GeV. 



In most phenomenological applications, the renormalization scheme of choice is MS and 
the preferred reference scale is 2 GeV. For this reason, we provide in Table [6] the values of 
the scheme dependent RCs Zs, Zp, Zt and Z q in the MS scheme at the scale of 2 GeV. The 
conversion from the RI-MOM scheme at fi = 1/a to MS at fi =2 GeV has been performed 
using renormalization group-improved perturbation theory at the N 2 LO for Z? and at the 
N 3 LO for Zs, Zp and Z q . For completeness, we also report in the table the results for the 
scale independent Zy and Za- Note, however, that the Ward identity determination of 
Zy, given in Tables |4] and |5j has a much better statistical accuracy. 

It should be remarked that in Table [6] we only quote the results obtained from RI- 
MOM renormalization conditions with the method Ml, since they are in general affected 
by smaller discretization effects and should hence bepreferably used in computations at 
fixed gauge coupling without continuum extrapolation^. The case of Fig. [TT]is an exception 
in this sense, because there an accidental cancellation (enhancement) of the cutoff effects 
coming from the bare pseudoscalar meson mass M ps and from Zp computed with the 
method M2 (Ml) leads to a similar scaling behaviour of the two sets of data corresponding 
to the M2 and Ml determinations of Zp. Moreover, with respect to the method Ml, the 
method M2 tends to give results that are affected by smaller statistical uncertainties, since 
the latter procedure does not involve the extrapolation to a 2 p 2 = of Eq. ( 136]) . 



4 Conclusions 

In this paper we have presented a non-perturbative determination of the RCs of bilinear 
quark operators for the tree-level Symanzik improved gauge action and the Nf = 2 twisted 
mass fermion action at maximal twist, which guarantees automatic (9(a)-improvement. 
The values of these constants do not depend on whether tm, OS or Wilson fermions are 
used. The RCs of the five bilinear quark operators and the quark field have been determined 
with the RI-MOM method. We have also obtained an independent estimate of Zy, based 
on a Ward identity, and have introduced a new method for Za and Zs/Zp, based on 
a combined use of standard tm and OS fermions. The main results for the RCs are 
collected in Tables HI [5] and O where they are also compared with the predictions of 

9 Anyway, the values of the RCs in the MS scheme at the scale of 2 GeV obtained with the method M2 
can easily be obtained for each a = a(/3) from their Ml-method counterparts and the RI-MOM results in 
Tables Hand [5] with the help of the relation Z^(2 GeV)| M 2/^r / (V a )|M2 = zf S (2 GeV)\ M1 /Zfi '(I/^Imi- 
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one-loop perturbation theory. We find that the differences between perturbative and non- 
perturbative determinations are large in the cases of Z$ and Zp. This result is of particular 
relevance to the lattice determinations of quark masses and the chiral condensate. 



Appendix: O(a)-improvement of RI-MOM renormali- 
zation constants 

In this Appendix we prove that: 

1. at maximal twist the RI-MOM "form factors" from which RCs are extracted [1] are 
automatically 0(a) improved at all p 2 values; 

2. at generic values of the twist angle, form factor improvement is guaranteed only at 
"large" p 2 values; i.e. where spontaneous chiral breaking effects can be neglected. 

In both cases the resulting renormalization constants will be improved. In the first case this 
is obvious because they are extracted from ratios of automatically improved form factors. 
In the second case the desired conclusion follows from the observation that RCs are short- 
distance quantities (insensitive to infrared effects), evaluated from large p 2 correlators. In 
the ensuing sections we give a detailed proof of these two statements. 

Improvement at maximal twist 

In the RI-MOM scheme the RC Z[ h of the bilinear operator 

O f T h = q f Tq h , (39) 

is defined by the condition 

lim Z^Zl h Tr[h{ h {p)P T ]\ 2 2 = 1, (40) 
iJ> h a -^n p =/* 

where /, h — 1,2 are flavour indices, and Pp is the projector onto the Dirac structure of 
the bilinear. The other quantities entering Eq. ( 14"U|) are defined by 



-i Tr^^(p)- 1 ] 

lim Z n 



1 (41) 



// pm7o1->o y 12 p 2 
with 

4 h (p) = S f q (p)- 1 Gl h (p)S^p)-\ (42) 
S>) = a 4 £e-^(gW(0)>, (43) 

X 

G{\p) = a 8 Yl e-^~^ (q f (x)O[ h (0)q h (y)) . (44) 



x.y 



Two useful observations are: (i) the "form factors" V q (p) = Tr[^ S^(p) x ] and Vp (p) = 
Pr] are H(4) invariant by construction. In the continuum limit they approach 
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quantities which are invariant under both 0(4) and parity, (ii) Z q carries no flavour index 
/ since, owing to its definition (j4ip . it depends quadratically on r q . 

The proof that our estimators of the RI-MOM RCs are 0(a)-improved at all quark 
masses and external momenta in maximally twisted QCD follows closely the argument 
of ref. [13] (subsequently refined in ref. [27]) for automatic improvement of correlation 
functions. It is based on the exact invariance of the maximally twisted QCD action under 
the transformation V x V d x — » — where \i q is the twisted quark mass. For 
completeness we recall the definition of the discrete transformations V (continuum parity) 
and T>d in the physical quark basis 

{Uo(x) — > Uo(x-p) , Uk{x) — > U\{x-p — ak) , k = 1, 2, 3 
qf{x) lo qt{x v ) (45) 
q f (x) -> g / (x P )7o 

f U^x) -^XjU-x-api) 

where x-p = (— x, xo) and (x is the unit vector in the /i-direction. Actually the proof is 
also valid in the more general "mixed action" setting, where Osterwalder-Seiler fermions 
are used as valence quarks. This is because the corresponding action is similarly invariant 
under V x T>d x (M q —> —M q ), where M q denotes the set of all quark (sea and valence) 
twisted mass parameters. 

The proof 

The proof can be given in three steps: 

• Due to H(4) symmetry, the lattice form factors V q (p) and V/ h (p) have no terms of 
the form p$ + p\ + p\ + p\ and (poPiP2P3) k \ with odd k. Thus they are invariant under 
parity transformation; i.e. they satisfy the relations 

V ff (p) = V q {p v ) , vl\p) = Vl\p v ) 1 p v = (- P ,p ). (47) 

• As has been proved in [27], the invariance of maximally twisted QCD under V x T>d x 
(iiq — > —fx q ) implies that 0(a J ) cut-off artefacts of lattice correlators arise, in their 
Symanzik expansion, from the insertion of Symanzik operators of parity (— 

• Subsequently, 0(a J ) terms with odd j must then be absent in the Symanzik expansion 
of parity-even form factors, like V q (p) and Vp h (p). 

The Symanzik expansion of the form factors above is straightforwardly obtained from that 
of the correlators in terms of which they are defined. 

Improvement at a generic value of the twist angle 

Since RCs are UV quantities, they should be independent of the QCD infrared structure 
and properties; in particular they should be unaffected by spontaneous chiral symmetry 
breaking. Thus, if we compute RCs at zero quark mass and at high momentum scale p 2 
such that chiral breaking effects are negligible, we can fully exploit the 7?.5-symmetry of 
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the continuum theory. This will be useful because in turn the lattice action is invariant 
under 1Z 5 x V d , where 

*.=n*f. m 

is a transformation belonging to the SU (Nf)LxS\J(Nf)n chiral group in the unitary case 
(or the graded SU(N s f ea + NJ al \NJ al ) L xSV(Nf ea + NJ al \NJ al ) R chiral group in the mixed 
action case). 

The proof 

The proof of the 0(a) improvement of RCs in lattice QCD at generic values of the twist 
angle, e.g. of the bare twist angle u = arctan(/i 9 /(m — M cr )), runs on lines similar to 
those of the previous section, with 1Z 5 replacing V . One can, in fact, argue as follows. 

• RCs constants can be extracted from the large-p 2 behaviour of the V q (p) and V[ h (p) 
form factors in the chiral limit. In this regime chiral breaking effects can be ignored and 
Vq(p) and V[ h (p) are invariant under 71$. 

• The invariance of the theory under 71$ x Dd implies that the ^-parity of an operator 
equals the parity of its naive dimension. Thus, 0(a J ) lattice artefacts of correlators or form 
factors arise in the Symanzik expansion from the insertion of operators of ^-parity equal 
to (-l)J'. 

• Therefore 0(a J ) terms with odd j cannot be present in the Symanzik expansion of IZ5- 
even form factors, like V q (p) and Vp (p). In fact, the continuum target theory is invariant 
under the chiral group (including IZ5), while 0(a J ') terms would be odd under TZ 5 . 

A final comment 

The 0(a) improvement of RI-MOM RCs at any value (including zero) of the twist angle uq 
hardly comes as a surprise because the RC-estimators computed at generic uo become uq- 
independent in the chiral limit, provided they are analytic in the complex mass parameter 
(m — M cr ) + ifi q , which is certainly the case at sufficiently large p 2 . As a consequence, a 
unique set of RI-MOM RCs, all free from 0(a^) (j odd) artifacts, exist for lattice QCD 
with Wilson quarks and a given pure gauge lattice action. These RCs can in principle 
be extracted from simulations at arbitrary (and not necessarily all equal) values of the 
twist angle. Computing the RI-MOM RCs at maximal twist (|a>o| = k/2), as we do in the 
present study, is just a simple and technically convenient choice. In our case this choice 
was a quite natural one, as we could then evaluate the necessary correlators on the Nf = 2 
gauge ensembles which were produced for studying large volume physics. 
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